plants <- read.table("ecol 563/jimsonweed.txt", header=T, sep='\t')
plants[1:12,]
 
# display layout of experiment
library(lattice)
dotplot(factor(pot)~lw.rat, groups=type, data=plants)
dotplot(factor(pot)~lw.rat, groups=type, data=plants, auto.key=T)
 
# fit randomized block design
mod1 <- lm(lw.rat~factor(pot)+type, data=plants)
anova(mod1)
summary(mod1)
 
# Model ignoring blocks: compare MSE to blocks model
mod0 <- lm(lw.rat~type, data=plants)
summary(mod0)
 
# test for a treatment by block interaction
mod2 <- lm(lw.rat~factor(pot)*type, data=plants)
anova(mod2)
 
# to predict the mean treatment we also need to know the block
predict(mod1, newdata=data.frame(type=c('G','N')))
 
# generate all unique combinations of pot and type
# type comes first to match the order of the data in data set
newdata <- expand.grid(type=levels(plants$type), pot=unique(plants$pot))
out.p <- predict(mod1, newdata=newdata)
out.p
cbind(newdata, out.p)
 
# treat blocks as random: using nlme package
library(nlme)
mod2.lme <- lme(lw.rat~type, random=~1|pot, data=plants)
 
# compare random blocks with fixed blocks model
anova(mod1)
anova(mod2.lme)
summary(mod2.lme)
coef(mod1)
 
# extract fixed effects estimates with fixef
fixef(mod2.lme)
# extract random effects predictions with ranef
ranef(mod2.lme)
# coef combines random effects and fixed effects
coef(mod2.lme)
 
# the predict function uses both fixed effects and random effects by default
predict(mod2.lme)[1:10]
 
# specifying level=0 gets predictions using just the fixed effects
predict(mod2.lme, level=0)[1:10]
 
# treat blocks as random: using lme4 package
library(lme4)
 
# to avoid function conflicts unload nlme from memory
detach(package:nlme)
 
# fit randomized block design
mod2.lmer <- lmer(lw.rat~type+(1|pot), data=plants)
 
# lmer does not return p-values
anova(mod2.lmer)
summary(mod2.lmer)
 
# The F-statistic is in the ANOVA table
anova(mod2.lmer)[,4]
 
# use a parametric bootstrap to obtain a p-value
# fit a model to the data without type as a predictor
mod1.lmer <- lmer(lw.rat~(1|pot), data=plants)
 
# start Monte Carlo simulation--1000 simulations
nrep <- 1000
#initialize storage vector
Fstat <- numeric(nrep)
 
# loop to carry out simulations
for(i in 1:nrep) {
 # simulate data from model in which type has no effect
 rmath <- unlist(simulate(mod1.lmer))
 # estimate type model to these data
 rmod <- lmer(rmath~(1|pot)+type, data=plants)
 # extract statistic
 Fstat[i] <- anova(rmod)[1,4]
}
 
max(Fstat)
# null distribution of F-statistic
hist(Fstat)
 
# p-value of actual F-statistic
1/1001
 
# Bayesian estimation using MCMC
out.mc <- mcmcsamp(mod2.lmer, n=10000)
# this creates an S4 object
slotNames(out.mc)
# convert to a data frame
ff <- as.data.frame(out.mc)
head(ff)
 
# posterior distribution of the treatment effect
hist(ff$typeN)
# mean of the posterior distribution
apply(ff, 2, mean)
fixef(mod2.lmer)
 
# 95% credible interval for treatment effect using percentile method
apply(ff, 2, function(x) quantile(x, c(.025,.975)))
 
# 95% higest probability density interval
HPDinterval(out.mc)
 
# obtain posterior distributions of the treatment means
mean.post <- data.frame(mean.g=ff[,1], mean.n=ff[,1]+ff[,2])
mean.post[1:10,]
 
# obtain 95% credible intervals for the individual treatment means
apply(mean.post, 2, function(x) quantile(x, c(.025,.975)))